SUBROUTINE vfunapprox_thetav_only(emaxin,betain)
  USE Commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'

  REAL(8), INTENT(IN)  :: emaxin(Ntheta_vot)
  REAL(8), INTENT(OUT) :: betain(MMc)
  REAL(8), ALLOCATABLE :: y(:), xvar(:,:), vfunout(:), ones(:), xpx(:,:), xpy(:), invxx(:,:), res_final(:), work(:)
  REAL(8)              :: Rsquared
  INTEGER              :: ii, j, itheta, info, MM_t, lwork
  CHARACTER :: uplo='U'

  ii = Ntheta_vot

  ALLOCATE(y(ii), ones(ii), vfunout(Ntheta_vot), res_final(ii))
  y = 0d0; res_final = 0d0

  ! In the following 1 is the constant
  MM_t = MMc
  ALLOCATE(xvar(ii,MM_t),xpx(MM_t,MM_t), xpy(MM_t), invxx(MM_t,MM_t))

  !We compute the coefficients on public consumption in V = f(sav)+alpha_1+alpha_2*q+alpha_3*q^2
  ones = 1d0
  DO itheta = 1, Ntheta_vot

     y(itheta) = emaxin(itheta)
     xvar(itheta,1:MM_t) = (/ ones(itheta), theta(itheta) /)
           
  END DO
  
  xpx = MATMUL(TRANSPOSE(xvar(:,:)),xvar(:,:))/SIZE(y)
  xpy = MATMUL(TRANSPOSE(xvar(:,:)),y(:))/SIZE(y)
     
  invxx = 0d0
  DO j=1, MM_t
     invxx(j,j)=1d0
  END DO

  LWORK = MM_t
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))  
  CALL DPOSV(uplo,MM_t,MM_t,xpx,MM_t,invxx,MM_t,info)
  IF (info /= 0) THEN
     DO j = 1, ii
        WRITE(*,'(a,1i3,7f12.6)') 'y,x',j,y(j),xvar(j,1:MM_t)
     ENDDO
     WRITE(*,*) 'MATRIX COULD NOT INVERT vfun__thetav_only 1', 'info', info
  END IF
  
  betain(:) = MATMUL(invxx,xpy)

  !Compare the actual value function with the approximation
  DO itheta = 1, Ntheta_vot

     vfunout(itheta) = betain(1) + betain(2)*theta(itheta)
     res_final(itheta) = emaxin(itheta) - vfunout(itheta)
     y(itheta) = emaxin(itheta)
           
  END DO

  Rsquared = (SUM((y-SUM(y)/SIZE(y))**2d0) - SUM(res_final**2d0)) / SUM((y-SUM(y)/SIZE(y))**2d0)

  IF (Rsquared < 0.75d0 .AND. flag_Tend == 0) THEN
!!$  IF (Rsquared < 0.75d0) THEN
!!$  IF (flag_here == 1) THEN
!!$  IF (flag_Tend == 1) THEN
     IF (myrank == 0) THEN
        write(*,*) 'Rsquared vfunapprox_Tend', Rsquared
        DO itheta = 1, Ntheta_vot
              
              write(*,'(a,2i4,5f20.10)') 'vfunapprox_thetav_only', itheta, emaxin(itheta), vfunout(itheta), betain(1) + betain(2)*theta(itheta)
              
        END DO
        write(*,*) 'betas', betain(:)
     END IF
     flag_here = 146
  END IF
  
  DEALLOCATE(y,xvar,res_final,vfunout)

END SUBROUTINE vfunapprox_thetav_only

